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Abstract  -  A  dynamic  analysis  of  the  Correlation  Integral  (Cm) 
of  the  Heart  Rate  Variability  signal  (HRV)  was  applied  to  50 
patients  with  Hypertrophic  Cardiomyopathy  (HCM).  A  group 
of  55  healthy  subjects  was  considered  as  a  control  group.  The 
Correlation  Integral  is  calculated  within  a  moving  time  window 
in  order  to  characterize  the  nonlinear  dynamical  behavior  of 
the  HRV  of  HCM  patients  that  cannot  be  described  by  classical 
correlation  dimension. 
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I.  Introduction 

The  analysis  of  Heart  Rate  Variability  (HRV)  has  been 
relevant  in  the  study  of  several  cardiovascular  phenomena. 
In  the  present  study  Hypertrophic  Cardiomyopathy  (HCM) 
patients  are  described  by  the  Heart  Rate  Variability  (HRV) 
using  complexity  analysis.  In  HCM  patients  there  is  an 
increased  risk  of  premature  death,  which  can  occur  with 
little  or  no  warning.  Sudden  cardiac  death  (SCD)  can  strike 
at  any  age. 

Several  studies  have  indicated  the  non-linear  nature  of 
heartbeat  dynamics  and  cardiac  function.  The  fluctuations  in 
time  between  beats  in  a  healthy  heart  may  be  due  to  an 
intrinsic  variability  of  the  autonomic  nervous  system  control 
of  non-linear  nature,  which  accounts  for  this  chaotic 
behavior,  while  illness  is  associated  to  an  increasing 
regularity  or  complexity  decrease  [1,2].  In  the  case  of 
hypertrophic  cardiomyopathy  patients,  some  attempts  to 
characterize  this  cardiac  illness  using  some  non-linear  tools 
like  Poincare  plots  and  entropies  have  been  made  [3,4].  One 
of  the  measures  of  the  time  series  complexity  is  the 
Correlation  Dimension  (Dc),  which  gives  a  statistical 
measure  of  the  self- similarity  of  the  geometry  of  a  set  of 
points  (attractor)  in  a  phase  space  [5]. 

Based  on  these  concepts,  this  work  presents  a  dynamic 
analysis  of  the  Correlation  Integral  (Cm)  applied  to  the  Heart 
Rate  Variability  in  order  to  characterize  HCM  patients.  In 
this  study,  Cm  was  calculated  on  a  moving  window,  which 
was  moved  along  the  RR  file  several  beats  successively  for 
all  calculations.  This  methodology  permits  to  analyze 
complexity  changes  throughout  time.  The  number  of  beats 
(N)  considered  was  common  to  both  groups  (HCM  and 
NRM),  in  order  to  obtain  comparable  results.  The 
correlation  dimension  based  on  Grassberger  and  Procaccia 
[6]  criteria  was  also  applied  to  these  signals  without 
considering  the  moving  window. 


Furthermore,  the  dynamic  analysis  of  the  Correlation 
Integral  (Cm)  was  applied  to  RR  series  of  a  patient  with 
cardiac  arrest  and  a  patient  with  aortic  stenosis. 


II.  Analyzed  data 

HRV  signals  were  obtained,  using  software  developed  by 
our  group  [7],  from  24-hour  ECG  recordings  sampled  at  256 
Hz  of  50  patients  with  HCM  of  the  NIC  database  (National 
Institute  of  Cardiology,  Warsaw).  All  patients  had 
QRS<  120ms,  no  conduction  abnormalities,  less  than 
100/24h  ectopic  beats  and  did  not  receive  medication  during 
recording.  A  group  of  55  normal  subjects  (NRM)  was 
considered  as  a  control  group 

In  order  to  emphasize  specific  aspects  of  the 
methodology,  signals  from  two  subjects  without  HCM  were 
also  considered:  CA,  patient  with  cardiac  arrest;  AS,  patient 
with  aortic  stenosis.  These  Holter  recordings  were  sampled 
at  128  Hz  and  RR  interval  extracted  using  Del  Mar  563 
Strata  Scan  software. 

Heart  Rate  Variability  depends  on  many  factors.  One  of 
them  is  patient  activity  during  recording.  To  achieve  the 
most  comparable  recording  conditions  and  the  best  signal 
quality,  we  have  selected  RR  data  series  segments  of  4  night 
hours  recorded  approximately  from  2  to  6  am,  while  patients 
were  sleeping. 


III.  Methodology 


From  a  defined  data  vector  (yi?  i=l,...,N),  a  m 
dimensional  phase  space  (embedding  dimension)  is 
constructed,  according  to  Takens  theorem  [8],  obtaining 


Xt  =  (yt,  yt+x,  yt+2x,...,yt+(m-i)x),  t=  1,  ...,N-(m-l)x  (1) 


where  x  is  the  lag,  expressed  as  a  number  of  beats. 

In  the  phase  space,  the  Correlation  Integral  which 
measures  the  number  of  points  Xj  that  are  correlated  with 
each  other  in  a  sphere  of  radius  r  around  the  points  Xi  can  be 
defined  by 


Cm(r)  =  : 


Nref  N-(m-l)T 
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-Xj) 


(2) 


Nref[N  — (m  — l)x  +  l]tf  H 

where  N  is  the  number  of  points,  and  0(z)  is  the  Heaviside 
function 
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0(z)  =  0  if  z  <  0 
0(z)  =  1  if  z  >  0 

The  distance  II  Xi  -  Xj  II  between  a  pair  of  points  in  the 
attractor  has  been  computed  as  the  euclidean  norm 


In  practice,  it  has  been  shown  that  it  is  sufficient  to  take 
randomly  only  10%  of  the  points  as  reference  points  (Nref)  to 
calculate  the  Correlation  Integral  [9]. 

For  each  reference  point,  the  distance  to  all  other  points  in 
the  attractor  is  calculated  and  Theiler’  s  correction  is  applied 
to  avoid  autocorrelation  effects,  excluding  the  four  nearest 
neighbor  points  in  time  [10]. 

When  log2(Cm(r))  is  plotted  versus  log2(r)  (Fig.  1),  the 
slope  of  the  resulting  straight  line,  determined  by  a  linear 
regression  at  low  r,  yields  the  correlation  dimension  Dc  as 
was  showed  by  Grassberger  and  Procaccia  [6]. 


Fig.  1.  Correlation  Integral  as  a  function  of  the  sphere  radius  (r) 
for  each  embedding  dimension. 

There  have  been  some  proposals  concerning  the  amount 
of  data  points  N  needed  to  estimate  Dc  [11,12].  This 
presents  limitations  when  the  Dc  is  calculated  on  a  moving 
window  of  reduced  RR  samples.  For  this  reason,  we  propose 
to  estimate  a  correlation  integral  (Cm)  in  a  moving  time 
window  instead  of  using  the  concept  of  Dc.  In  this  way, 
several  Cm(r)  are  computed  for  increasing  values  of  m,  and 
the  slopes  (Sc)  are  determined  from  a  scaling  region  of  the 
log-log  plot,  obtaining  a  sequence  of  Sc(m). 

To  obtain  the  evolution  of  the  slopes  Sc  for  increasing  m, 
a  graph  of  the  Sc(m)  is  plotted  and  an  exponential  model  of 
the  type 

Sc=Sc*(  1  -exp(-km))  (4) 

is  fitted.  The  values  of  Sc*  and  k  are  estimated  using  the 
Levenberg-Marquardt  method. 

The  saturation  value  of  the  curve  Sc  ,  is  the  estimated 
value  of  the  Cm  slope  at  high  m  of  the  signal  (Fig.  2). 


with  a  fitted  exponential  curve. 

This  exponential  model  permits  to  calculate  more 
accurately  the  correlation  integral  slope  at  high  m  values 
because  all  slopes  of  the  correlation  integrals,  for  all  the 
different  embedding  dimensions  m,  are  used  at  the  same 
time.  The  exponential  curve  fit  follows  the  asymptotic 
behavior  of  the  calculated  Sc(m).  Thus,  it  is  no  longer 
necessary  to  define  arbitrarily  where  the  Sc  curve  becomes 
sufficiently  constant  to  yield  a  saturation  value  for  Sc.  The 
error  calculated  as  the  mean  square  root  of  the  Sc  values  and 
the  fitted  curve  was  less  than  5%,  with  respect  to  the  Sc* 
value. 

Based  on  this  proposed  methodology  we  present  a 
dynamic  analysis  of  the  Correlation  Integral  (Cm)  using  a 
moving  window  on  the  RR  signal.  In  this  approach,  Cm  is 
calculated  for  a  moving  window  of  2,000  points,  which  is 
moved  100  samples  along  the  signal  successively  for  each 
new  calculation  of  the  Sc.  In  order  to  obtain  an  accurate 
measure  of  Sc  we  considered  embedding  dimensions  m=l, 
2,. .  .,20,  and  the  signal  as  belonging  to  a  system  with  x=5. 

In  order  to  emphasize  specific  aspects  of  the 
methodology,  the  dynamic  analysis  of  the  Correlation 
Integral  (Cm)  from  signals  of  CA  and  AS  patients  was  done 
on  the  total  night  period. 

The  number  of  beats  (N)  considered  were  common  to  both 
groups  (HCM  and  NRM)  of  subjects  considered,  N=10,000, 
in  order  to  obtain  comparable  results.  The  correlation 
dimension  based  on  Grassberger  and  Procaccia  [6]  criteria 
was  also  applied  to  these  signals  without  considering  the 
moving  window,  on  the  N= 10,000  beats  and  also  during  the 
total  4  night  hours. 

IV.  Results 

To  show  the  usefulness  of  window  correlation  integral  we 
present  the  results  obtained  from  RR  series  of  a  patient  who 
had  had  cardiac  arrest  (CA)  (Fig.  3).  As  can  be  seen  in  the 
tachogram  of  Fig.  3. a,  long  periods  of  high-level  arrhythmia 
are  broken  by  short  periods  of  extremely  low  variance 
behavior.  Except  for  a  slight  shift  in  time,  approximately  the 
length  of  the  time  window,  the  dependence  of  the  window 
correlation  integral  is  strongly  correlated  with  the  behavior 
of  the  tachogram  (Fig.  3.b).  Contrarily  to  what  one  may 
expect  from  the  shape  of  the  tachogram  for  this  case,  the 


3  of  4 


window  correlation  integral  slope  decreases  sharply  during 
periods  of  arrhythmia.  This  is  due  to  the  fact  that  in  3- 
dimensional  phase  space,  the  trajectory  for  these  periods  is 
extremely  simple  (Fig.  3.c).  On  the  other  hand,  the  phase 
space  trajectory  becomes  much  more  complex  during  the 
low  variance  HRV  periods  (Fig.  3.d). 


time  (h) 


Fig.  3.  (a)  RR  series  of  CA  patient,  (b)  Slopes  of  the  correlation  integral  at 
high  m  values,  obtained  with  windows  of  2,000  beats  moving  each  window 
100  beats,  (c)  Three-dimensional  phase  space  trajectory  of  200  beats  during 
the  high  variance  of  HRV.  (d)  Three-dimensional  phase  space  trajectory  of 
200  beats  during  the  low  variance  of  HRV. 

Fig.  4  presents  the  results  of  a  patient  with  aortic  stenosis 
(AS),  in  which  characteristic  spiking  in  the  HRV  occurs. 
This  spiking  is  unusually  regular  (Fig.  4. a).  As  a  result,  the 
windowed  correlation  integral  slope  is  low  and  shows  very 
little  variance  (Fig.  4.b).  The  3 -dimensional  phase  space 
trajectory  is  shown  in  Fig.  4.c. 


Fig.  4.  (a)  RR  series  of  AS  patient,  (b)  Slopes  of  the  correlation  integral  at 
high  m  values,  obtained  with  windows  of  2,000  beats  moving  each  window 
100  beats,  (c)  Three-dimensional  phase  space  trajectory  of  2,000  beats. 

Fig.  5  shows  the  evolution  through  the  time  of  the  slopes 
(Sc)  at  high  embedding  dimension  m  of  the  correlation 
integral  (Cm)  of  a  RR  signal  behavior,  of  a  HCM  patient. 


Fig.  5.  (a)  RR  series  of  a  HCM  patient,  (b)  Evolution  of  the  slopes  (Sc)  of 
the  windowed  correlation  integral  at  high  embedding  dimensions. 

(c)  Three-dimensional  phase  space  trajectory  of  2,000  beats. 


Table  I  presents  the  results  of  applying  the  proposed 
methodology  to  the  HCM  and  NRM  groups.  The  mean 
value  and  the  standard  deviation  (sd)  of  Sc*  obtained  in  all 
studied  windows  characterized  HCM  group  (mean  Sc*: 
10.0+0.955;  sd  Sc*:  1.589+0.624)  with  a  significant  level 
p=0.041  and  p<0.0005,  respectively,  when  both  groups  were 
compared.  The  minimum  value  of  Sc*  could  differentiate  the 
groups  (HCM:  6.603+1.755;  NRM:  7.472+1.594)  with  a 
p=0.006.  No  statistical  significant  differences  were  found 
considering  the  maximum  values  of  Sc*,  however  the 
differences  between  the  maximum  and  minimum  value 
(HCM:  4.876+1.409;  NRM:  5.949+1.807)  presented  a 
significant  level  p=0.001.  The  sum  of  the  consecutive 
absolutes  differences  of  the  Sc*  values  were  also  statistically 
significant  with  p=0.026. 

Table  i 


SLOPES  (Sc)  AT  HIGH  EMBEDDING  DIMENSION  OF  THE 
WINDOWED  CORRELATION  INTEGRAL  Cm 


mean+sd 

groups 

NRM 

HCM 

p-value 

mean  Sc 

10.4+0.786 

10.0+0.955 

0.041 

sd  Sc* 

1.186+0.456 

1.589+0.624 

<0.0005 

min  Sc* 

7.472+1.594 

6.603+1.755 

0.006 

max  Sc* 

12.3+0.725 

12.6+0.716 

n.s. 

maxSc*-minSc* 

4.876+1.409 

5.949+1.807 

0.001 

IISc*(i)-Sc*(i+l)l 

0.389+0.067 

0.426+0.089 

0.026 

n.s.  no  significant  statistical  level 


It  can  be  seen  in  Table  II  that  Correlation  Dimension  (Dc) 
calculated  from  a  non-windowed  correlation  integral  and 
considering  the  complete  night  period  of  the  RR  series  could 
not  characterize  the  system  complexity  of  HCM  patients. 
Similarly,  when  using  Dc  to  analyze  10,000  samples  of  the 
RR  series  (Table  III),  HCM  and  NRM  subjects  were  not 
differentiated.  However,  both  groups  could  be  characterized, 
in  the  time  domain,  by  the  mean  value  of  the  RR  series  and 
the  standard  deviation  (Tables  II  and  III).  Moreover,  no 
statistical  differences  were  found  comparing  Dc  from  the 
complete  night  period  and  10,000  samples  of  the  RR  series. 
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Table  ii 


CORRELATION  DIMENSION  Dc  OF  THE  RR  SERIES  DURING 
4  HOUR  NIGHT  PERIOD 


mean+sd 

groups 

NRM 

HCM 

p-value 

Dc 

9.721+1.346 

9.249+1.896 

n.s. 

mean  RR 

960.2+132.3 

1028.9+154.3 

0.028 

sd  RR 

100.9+53.4 

114.8+41.8 

0.023 

n.s.  no  significant  statistical  level 

Table  hi 

CORRELATION  DIMENSION  Dc  OF  THE  RR  SERIES 
10,000  SAMPLES 

mean+sd 

groups 

NRM 

HCM 

p-value 

Dc 

9.503+1.196 

9.013+1.781 

n.s. 

mean  RR 

958.3+131.1 

1027+155.7 

0.016 

sd  RR 

101.0+45.9 

116.1+40.6 

0.026 

n.s.  no  significant  statistical  level 


V.  Discussion  and  conclusions 

A  methodology  to  obtain  dynamic  Correlation  Integral  of 
the  HRV  signal  has  been  presented.  Using  this 
methodology,  50  patients  with  Hypertrophic 
Cardiomyopathy  and  55  subjects  considered  as  a  control 
group  were  analyzed. 

With  respect  to  the  methodology,  it  is  important  to 
emphasize  that  the  exponential  fit  to  obtain  the  final  value  of 
the  Correlation  Integral  slope,  named  here  Sc*,  seems  to  be  a 
better  approach  to  obtain  the  correct  value,  since  it  uses  all 
the  points  of  the  Sc(m)  curve  and  defines  the  correlation 
integral  slope  as  an  asymptotic  value  attained  with 
increasing  embedding.  The  curve  may  be  fit  with  a  relative 
low  error  -  less  than  5%  of  the  group  of  HCM  patients 
studied  here.  With  regard  to  the  number  of  beats  to  analyze, 
it  seems  that  10,000  beats  give  sufficient  information  to 
calculate  the  evolution  of  Sc*  consistently,  if  the  night  hours 
are  chosen  correctly. 

As  many  biological  phenomena  are  complex  the  approach 
of  moving  window  correlation  integral  Cm  could  be  an 
interesting  way  to  analyze  such  systems.  The  value  of  Sc* 
for  each  window  usually  fluctuates  around  an  average  value 
but,  as  may  be  clearly  seen  for  the  arrhythmia  case  CA,  may 
give  valuable  information  about  the  changes  in  the  dynamics 
of  the  system.  In  particular,  window  correlation  integral  may 
be  used  to  detect  relatively  short  episodes  embedded  in  the 
time  series.  In  the  case  of  HCM  patients,  changes  in 
complexity  evaluated  by  the  windowed  correlation  integral 
seems  to  be  able  to  characterize  this  group  of  patients. 
Contrarily,  the  results  obtained  from  HCM  patients  applying 
correlation  dimension  Dc  without  considering  a  moving 
window,  shows  that  Dc  does  not  characterize  heart  rate 
variability  in  HCM  patients. 
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